﻿///////////////////////////////////////////////////////////////////////////////
//
//  This file is part of MathLib.NET.
//
//  This library is free software; you can redistribute it and/or
//  modify it under the terms of the GNU Lesser General Public
//  License as published by the Free Software Foundation; either
//  version 2.1 of the License, or (at your option) any later version.
//  
//  This library is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//  Lesser General Public License for more details.
//  
//  You should have received a copy of the GNU Lesser General Public
//  License along with this library;  If not, see 
//  <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////////

using System;
using ILNumerics;
using ILNumerics.BuiltInFunctions;

// Math Operations
// ===============
//
// divide - Divide two objects 

namespace MathLib
{
    /// <summary>
    /// Mathematical constants.
    /// </summary>
    public partial class MLMath
    {
        //
        // --- Divide two objects ---
        //

        /// <summary>
        /// Divide two objects.
        /// </summary>
        /// <param name="op1">First operand (numerator)</param>
        /// <param name="op2">Second operand (denominator)</param>
        /// <returns>Division of objects (numerator/denominator)</returns>
        public static ILArray<complex> divide(ILArray<complex> op1, ILBaseArray op2)
        {
            return ILMath.divide(op1, ILMath.tocomplex(op2));
        }

        /// <summary>
        /// Divide two objects.
        /// </summary>
        /// <param name="op1">First operand (numerator)</param>
        /// <param name="op2">Second operand (denominator)</param>
        /// <returns>Division of objects (numerator/denominator)</returns>
        public static ILArray<fcomplex> divide(ILArray<fcomplex> op1, ILBaseArray op2)
        {
            return ILMath.divide(op1, ILMath.tofcomplex(op2));
        }

        /// <summary>
        /// Divide two objects.
        /// </summary>
        /// <param name="op1">First operand (numerator)</param>
        /// <param name="op2">Second operand (denominator)</param>
        /// <returns>Division of objects (numerator/denominator)</returns>
        public static ILArray<double> divide(ILArray<double> op1, ILBaseArray op2)
        {
            return ILMath.divide(op1, ILMath.todouble(op2));
        }

        /// <summary>
        /// Divide two objects.
        /// </summary>
        /// <param name="op1">First operand (numerator)</param>
        /// <param name="op2">Second operand (denominator)</param>
        /// <returns>Division of objects (numerator/denominator)</returns>
        public static ILArray<float> divide(ILArray<float> op1, ILBaseArray op2)
        {
            return ILMath.divide(op1, ILMath.tosingle(op2));
        }

        /// <summary>
        /// Divide two objects.
        /// </summary>
        /// <param name="op1">First operand (numerator)</param>
        /// <param name="op2">Second operand (denominator)</param>
        /// <returns>Division of objects (numerator/denominator)</returns>
        public static ILArray<long> divide(ILArray<long> op1, ILBaseArray op2)
        {
            return ILMath.divide(op1, ILMath.toint64(op2));
        }

        /// <summary>
        /// Divide two objects.
        /// </summary>
        /// <param name="op1">First operand (numerator)</param>
        /// <param name="op2">Second operand (denominator)</param>
        /// <returns>Division of objects (numerator/denominator)</returns>
        public static ILArray<ulong> divide(ILArray<ulong> op1, ILBaseArray op2)
        {
            return ILMath.divide(op1, ILMath.touint64(op2));
        }

        /// <summary>
        /// Divide two objects.
        /// </summary>
        /// <param name="op1">First operand (numerator)</param>
        /// <param name="op2">Second operand (denominator)</param>
        /// <returns>Division of objects (numerator/denominator)</returns>
        public static ILArray<int> divide(ILArray<int> op1, ILBaseArray op2)
        {
            return ILMath.divide(op1, ILMath.toint32(op2));
        }

        /// <summary>
        /// Divide two objects.
        /// </summary>
        /// <param name="op1">First operand (numerator)</param>
        /// <param name="op2">Second operand (denominator)</param>
        /// <returns>Division of objects (numerator/denominator)</returns>
        public static ILArray<uint> divide(ILArray<uint> op1, ILBaseArray op2)
        {
            return ILMath.divide(op1, ILMath.touint32(op2));
        }

        /// <summary>
        /// Divide two objects.
        /// </summary>
        /// <param name="op1">First operand (numerator)</param>
        /// <param name="op2">Second operand (denominator)</param>
        /// <returns>Division of objects (numerator/denominator)</returns>
        public static ILArray<short> divide(ILArray<short> op1, ILBaseArray op2)
        {
            return ILMath.divide(op1, ILMath.toint16(op2));
        }

        /// <summary>
        /// Divide two objects.
        /// </summary>
        /// <param name="op1">First operand (numerator)</param>
        /// <param name="op2">Second operand (denominator)</param>
        /// <returns>Division of objects (numerator/denominator)</returns>
        public static ILArray<ushort> divide(ILArray<ushort> op1, ILBaseArray op2)
        {
            return ILMath.divide(op1, ILMath.touint16(op2));
        }

        /// <summary>
        /// Divide two objects.
        /// </summary>
        /// <param name="op1">First operand (numerator)</param>
        /// <param name="op2">Second operand (denominator)</param>
        /// <returns>Division of objects (numerator/denominator)</returns>
        public static ILArray<double> divide(ILArray<byte> op1, ILBaseArray op2)
        {
            return ILMath.divide(ILMath.todouble(op1), ILMath.todouble(op2));
        }

        /// <summary>
        /// Divide two objects.
        /// </summary>
        /// <param name="op1">First operand (numerator)</param>
        /// <param name="op2">Second operand (denominator)</param>
        /// <returns>Division of objects (numerator/denominator)</returns>
        public static ILArray<double> divide(ILArray<char> op1, ILBaseArray op2)
        {
            return ILMath.divide(ILMath.todouble(op1), ILMath.todouble(op2));
        }


        //
        // --- Matrix power ---
        //

        /// <summary>
        /// Matrix power
        /// </summary>
        /// <param name="op1">First operand</param>
        /// <param name="op2">Second operand</param>
        /// <returns>First operand raised to power of second operand</returns>
        public static ILArray<complex> mpower(ILBaseArray op1, ILBaseArray op2)
        {
            if (op1.Dimensions[0] != op1.Dimensions[1] || op2.Dimensions[0] != op2.Dimensions[1])
                throw new MathException("mpower(): Matrices must be square.");

            ILArray<complex> result;

            // Convert inputs to complex arrays
            ILArray<complex> op1c = ILMath.tocomplex(op1);
            ILArray<complex> op2c = ILMath.tocomplex(op2);

            if (op2.Dimensions[0] == 1)
            {
                // Raise matrix to scalar power
                result = new ILArray<complex>(op1.Dimensions.ToIntArray());

                // Copy input matrix base
                MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix m = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(op1.Dimensions[0], op1.Dimensions[1]);
                for (int r = 0; r < op1.Dimensions[0]; r++)
                    for (int c = 0; c < op1.Dimensions[1]; c++)
                        m[r, c] = new System.Numerics.Complex(
                            op1c.InternalArray4Experts[c * op1.Dimensions[0] + r].real,
                            op1c.InternalArray4Experts[c * op1.Dimensions[0] + r].imag);

                // Copy input scalar power
                System.Numerics.Complex n = new System.Numerics.Complex(
                    op2c.InternalArray4Experts[0].real,
                    op2c.InternalArray4Experts[0].imag);

                // Compute power
                MathNet.Numerics.LinearAlgebra.Complex.Factorization.DenseEvd evd = new MathNet.Numerics.LinearAlgebra.Complex.Factorization.DenseEvd(m);
                MathNet.Numerics.LinearAlgebra.Generic.Matrix<System.Numerics.Complex> evects = evd.EigenVectors();
                MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix P = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(evects.RowCount, evects.ColumnCount);
                for (int r = 0; r < evects.RowCount; r++)
                    for (int c = 0; c < evects.ColumnCount; c++)
                        P[r, c] = evects[r, c];
                MathNet.Numerics.LinearAlgebra.Generic.Vector<System.Numerics.Complex> evals = evd.EigenValues();
                MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix Dn = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(evals.Count, evals.Count);
                for (int i = 0; i < evals.Count; i++)
                    Dn[i, i] = System.Numerics.Complex.Pow(evals[i], n);
                MathNet.Numerics.LinearAlgebra.Generic.Matrix<System.Numerics.Complex> An = P * Dn * P.Inverse();

                // Copy to output
                for (int r = 0; r < op1.Dimensions[0]; r++)
                    for (int c = 0; c < op1.Dimensions[1]; c++)
                        result.InternalArray4Experts[c * op1.Dimensions[0] + r] = new complex(
                            An[r, c].Real, An[r, c].Imaginary);
            }
            else if (op1.Dimensions[0] == 1)
            {
                // Raise scalar to matrix power
                result = new ILArray<complex>(op2.Dimensions.ToIntArray());

                // Copy input scalar matrix
                System.Numerics.Complex n = new System.Numerics.Complex(
                    op1c.InternalArray4Experts[0].real,
                    op1c.InternalArray4Experts[0].imag);

                // Copy input matrix power
                MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix m = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(op2.Dimensions[0], op2.Dimensions[1]);
                for (int r = 0; r < op2.Dimensions[0]; r++)
                    for (int c = 0; c < op2.Dimensions[1]; c++)
                        m[r, c] = new System.Numerics.Complex(
                            op2c.InternalArray4Experts[c * op2.Dimensions[0] + r].real,
                            op2c.InternalArray4Experts[c * op2.Dimensions[0] + r].imag);

                // Compute power
                MathNet.Numerics.LinearAlgebra.Complex.Factorization.DenseEvd evd = new MathNet.Numerics.LinearAlgebra.Complex.Factorization.DenseEvd(m);
                MathNet.Numerics.LinearAlgebra.Generic.Matrix<System.Numerics.Complex> evects = evd.EigenVectors();
                MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix P = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(evects.RowCount, evects.ColumnCount);
                for (int r = 0; r < evects.RowCount; r++)
                    for (int c = 0; c < evects.ColumnCount; c++)
                        P[r, c] = evects[r, c];
                MathNet.Numerics.LinearAlgebra.Generic.Vector<System.Numerics.Complex> evals = evd.EigenValues();
                MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix Dn = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(evals.Count, evals.Count);
                for (int i = 0; i < evals.Count; i++)
                    Dn[i, i] = System.Numerics.Complex.Pow(n, evals[i]);
                MathNet.Numerics.LinearAlgebra.Generic.Matrix<System.Numerics.Complex> An = P * Dn * P.Inverse();

                // Copy to output
                for (int r = 0; r < op2.Dimensions[0]; r++)
                    for (int c = 0; c < op2.Dimensions[1]; c++)
                        result.InternalArray4Experts[c * op2.Dimensions[0] + r] = new complex(
                            An[r, c].Real, An[r, c].Imaginary);
            }
            else
            {
                // Raise matrix to matrix power is not allowed
                throw new MathException("mpower(): Cannot raise matrix to power of matrix.");
            }

            return result;
        }
    }
}
